Quantum percolation in granular metals 
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Theory of quantum corrections to conductivity of granular metal films is developed for the realistic 
case of large randomly distributed tunnel conductances. Quantum fluctuations of intergrain voltages 
(at energies E much below bare charging energy scale Ec) suppress the mean conductance g(E) 
much stronger than its standard deviation cr(E). At sufficiently low energies E* any distribution 
becomes broad, with cr(E t ) ~ g(E*), leading to strong local fluctuations of the tunneling density 
of states. Percolative nature of metal-insulator transition is established by combination of analytic 
and numerical analysis of the matrix renormalization group equations. 

PACS numbers: 71.30.+h, 64.60.Ak, 73.23.Hk 
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Introduction. — Low-temperature electron transport in 
granular metals was intensively studied during last 
years 0i S 0- It was shown that in the tempera- 
ture range T > gd (where g 3> 1 is characteristic value of 
dimensionless intergrain conductance in units of e 2 /2irh 
and 5 is the intragrain level spacing), quantum correc- 
tions to conductivity originate mainly from local fluc- 
tuations of voltages between neighboring grains. This 
effect can be treated within Ambegaokar-Eckern-Schoen 
model and leads to logarithmic temperature depen- 
dence of the effective conductance : 



9{T) = go 



4. g E c 
- hi — — — . 

z T 



(1) 



where Ec S> $ is the charging energy of an individual 
grain, go is the bare tunneling conductance of intergrain 
junctions (identical for all junctions), and z is the coor- 
dination number of the lattice [|| . The result Q is valid 
as long as the renormalized conductance g(T) is large, 
i.e. down to temperatures T\ = goEce~ Z9 °/ 4 . It was 
argued 0, that transition from metal to insulator be- 
havior (MIT, for brevity) occurs at T ~ T\ as long as 
T\ > gS. The same conclusion for the two-dimensional 
(2D) case was reached Q via instanton analysis. 

Although the above results may well be applied to ar- 
tificial 2D arrays of well-defined tunnel junctions, tunnel 
conductances gtj are random in granular metals. In this 
Letter we investigate the role of gij randomness for en- 
ergy (temperature) dependent properties of thin granular 
films, such as macroscopic conductance g c ^(T) and the 
local tunneling density of states (LTDoS) Vi(E). Quan- 
tum fluctuations lead to suppression of g^ described by 
the one-loop renormalization group (RG) equation: 



d 9ij _ 0n p 
^ ~ ^gij^j, 



(2) 



where t(E) = \n(g Ec/E) is the auxiliary RG "time", 
g being some mean bare conductance, and Rij is the 
resistance of the network between the points i and j. 
Physically, renormalization of g^ is due to fluctuations 




FIG. 1: Conducting bonds (with gij > 0) at different values 
of t. Results of numerical simulation of Eq. JSJ on the lattice 
20x20 with P (g) = {2g 6(g)/-K)/(g 2 +gl). a) t = 0.64g with 
the fraction of conducting bonds N con d = 0.95; b) t = 1.08go 
with iV cond = 0.55. 



of voltage between the grains i and j, which are governed 
by the corresponding resistance Rij . Eq. is a straight- 
forward generalization of the RG equation for a regular 
array 1] with g^ = g and = 2/gz, whose solution 
is given by Eq. QJ. The system of RG equations (J2J is 
nontrivial since R^ is a complicated nonlocal function of 
all individual conductances g^\. 

In a regular system, Eq. |J5| drives all conductances 
to zero simultaneously at t = t c = zg$/A, marking the 
point of the MIT. We will show that in a random system 
renormalized conductances of the junctions collapse to 
zero neither all at once, nor one by one, but in groups. 
These groups enclose clusters, consisting of one or sev- 
eral sites, which become disconnected from the rest of the 
network after the collapse (see Fig. QJ. As a result, the 
MIT in a natural granular system is a percolative tran- 
sition: it takes place, when enough clusters have become 
disconnected so that the percolation via still conducting 
links is destroyed. 

The above picture of conductances, eventually collaps- 
ing to zero, follows from the one- loop RG equation J2J. 
The one-loop approximation breaks down at g ~ 1; at 
lower energies the conductance decays exponentially with 
the RG "time" t(E). Therefore, Eq. J3J can adequately 
describe only the initial stage of the MIT. Nevertheless, 
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there exists a region near the transition where the per- 
colative cluster contains good conductances with g > 1 
so that Eq. is still applicable. 

We start from the case of relatively narrow original 
distribution Po(g) characterized by the mean value g 
and the standard deviation oo ^ ffo> an d show that the 
renormalizcd distribution P(g) broadens. In particular, 
for the square lattice (z = 4): 



a(E) cr 



9o/g 



9(E) g Q y/2]n(g /g)]n\n(g /gy 



(3) 



where g = g(E) = g Q — ln(g Ec/E). Eq. © is a large- 
ln(g /g) asymptotics of a more general expression (see 
Eq. ifTTJl below). It is valid as long as cr(E) <C g(E), i.e. 
above E* = T* = g E c e^ + ' To > T t , where T x marks 
MIT in an ideal array with cro <C 1. Thus transition 
from metal into insulator in a granular array is intrinsi- 
cally inhomogeneous. The vicinity of this transition at 
max(_E, T) < is difficult for exact analytical treatment 
as the width of distribution P(g\E) becomes of order of 
its mean value. In this region we employ the effective- 
medium approximation (EMA) and numerical solution of 
the RG equations J2J), and demonstrate that MIT is of 
percolative nature. 

Strong self-developed inhomogeneity of a granular ar- 
ray can be probed by scanning tunneling measurement of 
the LTDoS modified by the Coulomb zero-bias anomaly 
(ZBA) 0|1,@- ZBA modification Z{E) = v(E)jv a of 
the average LTDoS in a regular array was considered in 
Refs. P, y] and found to become very large before ap- 
proaching MIT. Here we analyze spatial fluctuations of 
the ZBA suppression factor Zi{E). For an originally nar- 
row distribution Po(g), the log- normal distribution of the 
ZBA factors is found, with std[ln Z t (E)\ w a(E)/g(E). 
Thus we predict order-of-unity local fluctuations of LT- 
DoS at ma,x(E,T) < T*. Spatial correlation length 
f (E) of these fluctuations was found to grow moderately 
with E decrease in the case of weak original disorder: 
((E) Ri y/\n\g Q /g(E)], reaching y/ln(g /ao) at the bor- 
der of strong inhomogeneity E ~ T*. For the region in 
the vicinity of MIT, where relative fluctuations are large, 
we present numerical analysis of LTDoS fluctuations. Be- 
low we provide brief derivation of our results. 

Narrow distribution. — If the standard deviation a of 
the distribution is much smaller than the mean g, the lat- 
ter follows the homogeneous solution 0J: g(t) = ~g Q — t, 
while evolution of Sgij = g^ — g can be described pertur- 
batively. Resistance can be written as Rij = Gu + Gjj — 
2Gij, where Gij = is the Green function of the dif- 
fusion operator on the network defined by the matrix 



elements A*. 



and Aij = 



-9ij 



L Dy 
10]. 



Using the 



standard perturbative series Gij 
we find 



Gij — Gih5AhiG 



6Gi 



(ki) 



Sgki(Gik — Gu)(Gjk — Gji) 



(4) 



where in momentum representation G(p) — [2'ge(p)]~ 1 ; 
for the square lattice, e(p) = 2 — cosp x — cosp y . 

To proceed further we choose a vector representation 
for conductances gf when each edge is characterized by 
the lattice site i it goes from and direction a which 
can be either horizontal (+x) or vertical (+y). Using 
Eq. introducing a new time variable s = hi[g /g(t)] — 
— ln[l — t/~g ] and passing to Fourier representation we 
get a linear evolution: 



dSg a (p) 
ds 



-Wl a p(p)Sg (p), 



governed by the 2x2 time-independent matrix 
3R(p) = 1 - Vi (p) - V 2 (p)(e la ^+ + e- M p<7_ 



(5) 



(6) 



where a p = (p x —p y )/2 and a± = (a\ T 02)/2, a k being 
the Pauli matrices. The functions Vi^(p) are given by 

Vi = 2 



(dq) (l-COBfa)(l-COBfe,-fa)) i (?a) 



^(q)e(p- q) 

— — COS q x ) (cos ?y 



, , s (cos — cos q x 1 1 cos '-if- — cos q v ) 
^ = 2/ (dqy- 2 . ■ q ;'\ 2 qy \ (7b) 



£ (q + p/2) £ (q-p/2) 

where the integral with (dq) = d 2 q/(2n) 2 runs over the 
Brillouin zone. At small p they have a nonanalytic be- 
havior: Vi{p) = l-l/ir - (p 2 /87r)ln(l/p) + ... and 
P 2 (p) = l/7r-(p 2 /87r) ln(l/p) + ... 

The eigenvalues of the matrix © form two branches: 



A±(p) = l-Pi(p)±7> 2 (p) ) 



(8) 



the eigenfunctions being (e Mp , =Fl) T /v / 2- The spectral 
branch A+(p) is gapped whereas the branch A_(p) be- 
comes gapless in the long wave-length limit: A_(p — > 
0) « (p 2 /An) ln(l/p). Once the spectral properties of the 
matrix 971 arc known one can express Sg's at time s via 
their initial values at s = 0: 



Sg a {r,s) = ^2K af3 (r-r', s)<%j(r',0). 



The Fourier-transformed kernel is given by 



(9) 



K(p, s) = Kx{p, s) + K 2 (p, s)(e ia »a+ + e^a^), (10) 

where K h2 (p,s) = ( e - x -^ s ± e~ x + (p >)/2. 

Equation JjJJ allows to find the evolution of the single- 
site distribution function P(g). A more convenient 
quantity is the characteristic function x(A) defined as 
e -x(A) = j p(Sg)e lXSg d5g. Assuming that at s = dif- 
ferent conductances are uncorrelated we find for x(A): 
x(A;s) = Er(x[^i(r,5)A;0] + x[K 2 (r,s)X;0]). Thus, 
the variance of the distribution will decay as 



a 2 (s) 



(dp) 



-2\-(p)s _|_ e -2A+(p) S 



(11) 
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FIG. 2: Evolution of a(s)/ao, Eq. 1111 . (dashed line), the re- 
sult of numerical simulation of Eq. with ffo/g = 0.1 (solid 
line), and the EMA prediction e~ a ^ 2 (dotted line). Inset: g(t) 
and er(t) near the MIT. 



For s <C 1, cx{s)/(jq = 1 — s/2 + . . . This series is 
to be compared with decay of the mean conductance: 
g(s)/~g Q — e~ s = 1 — s + . . . . Thus, even at the initial 
stage of the evolution the width of the distribution decays 
slower than its average. In the case s ^> 1 the integral 
(|llfl is dominated by the soft mode A_(p) at p — > lead- 
ing to cr 2 (s)/cr 2 w 1/ (2s Ins) and hence to Eq. (J3J). Prac- 
tically, the applicability of this asymptotics is limited to 
very large s. For intermediate values of s one has to 
employ full Eq. (|llfl with numerical integration over the 
Brillouin zone. The obtained function <j(s)/(Jq together 
with the prediction of the EMA and results of numerical 
simulation of Eq. J2J [for a typical realization of disorder 
on the 20 x 20 lattice with periodically boundary con- 
ditions and Po(g) — (2g 0(g)/n)/(g 2 + g$)] is shown in 
Fig.H 

Apart from broadening the single-site distribution 
P(g), the RG flow J2J produces correlations between 8g 
at different links: C a p(r] s\, s 2 ) = (5g a (r, S\) 8gp{0, 52))- 
The Fourier transform of the correlation function reads: 



C(p; si, s 2 ) = a 2 K(p; s i + s 2)- 



(12) 



At the initial stage of evolution, at s = (si + 
s 2)/2 ^ f, correlations are short-ranged. At the 
later stage, s > 1, correlations with large correlation 
length £(s) = \J (4/7r)slns develop: C a p(r, s%, S2) — 
CT 2 (s)exp[-r 2 /£ 2 (s)]- 

Spatial fluctuations of gij lead to fluctuations of the 
LTDoS Vi{E) = Zi{E)v Q . The ZBA suppression factor 
Zi(E) for granular media at E > gS can be found ac- 
cording to simple "environmental theory" : 



\nZi(E) 



Ri{t')d£ 



(13) 
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FIG. 3: Histograms for the distribution of the ZBA factors 
Z(E) near the percolation threshold t c = 0.99 ~g obtained 
numerically for <To/g = 0.32: t — 0.8 t c (solid line); t = 0.9 t c 
(dashed line); t = 0.95 t c (dotted line). 



region of the array at the energy scale E = ~g Vl Ece~ t . The 
same result follows from the analysis provided in 0, 0] . 
It is important to note that in a homogeneously disor- 
dered metal the short-length cutoff in the integral that 
determines the effective resistance R{E) is given by the 
diffusion length y/HD/E, whereas in the present case it 
is just the grain size. The long-scale cutoff for logarith- 
mic divergency of Ri(E) in 2D is L(E) — e 2 g c s/E (in 
the absence of external screening). Thus one can write 
Ri{E) = G-° g , where the otherwise divergent Ga is reg- 
ularized by the finite length L(E). Local fluctuations of 
In Zi(E) are determined by a much smaller region of the 
size £ around the site i so that the object 5Ri = SGu is 
already free of infra-red divergency and is independent 
on the details of screening. Employing Eq. Q we obtain 



ShxZi(t) = -2 



5g a (r k ,t') 



Uvi-Tk), (14) 



where Q a {v) is specified by its Fourier transform: 



Q a (p) 



- p -iPc./2 



rj„\ cos(p a /2) - cosg a 
(dq) 2e(q + p/2Mq-p/2)- (15) 



Averaging 6 In Zi(t) with the help of Eq. Ijl2|l . and chang- 
ing integration variable from t to s we get for the variance 
of the ZBA exponent: 



where Ri (t) is the resistance between the site i and the far 



([SlnZ^t)} 2 ) =4=1 / / d Sl ds 2 e Sl+S2 
So Jo Jo 

x J(dp)K afj (p;s 1 +s 2 )Q a (p)Q/3(-p). (16) 

At the initial stage of the evolution, for s < 1, 

([5lnZi(t)} 2 ) w (0.26cr /g ) 2s2 - In thc rc g ion of un- 
developed correlations, for s > 1, Eq. fro| is domi- 
nated by small momenta. In this limit Eq. I|15|) can 
be estimated as Q a (p — > 0) = (l/47r) ln(l/p), yielding 
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{[5\nZi(t)} 2 ) « (cr 2 /8n 2 g 2 )(e 2s \ns/s). Rewriting this 
result as ([8 In Zi(t)] 2 ) w (lns/27r) 2 cr 2 (s)/g 2 (s) we see 
that fluctuations of the ZBA factors Zi become of the 
order of unity simultaneously with the renormalized ra- 
tio a/g. 

The results of numerical simulation for a model dis- 
tribution P(g) = cxp[— (lng/g ) 2 /2a 2 ]/(V2Trcrig) with 
the moderately small variance cr 2 . = g 2 \e 2<Jl — e ai ] = 
(0.32 g ) 2 on the lattice 32 x 32 are shown in Fig. 
where we present the distribution of the local values of 
Zi{E) at three values of the RG "time" t. Upon lower- 
ing the energy scale and approaching the MIT transition 
at E ~ T c — g Ece~ tc with t c = 0.99g , we observe a 
growing relative width of Z distribution, with the zevo-Z 
peak developing near the percolation threshold, due to 
considerable weight of disconnected clusters. 

Effective medium approximation (EMA). — In this ap- 
proximation one takes into account only the simplest - 
local - correlations between <?y and i2y , while all distant 
conductances are replaced with a homogeneous medium 
with effective conductance g e $ (see, e.g., [HI). Spatial 
correlation are neglected within EMA, and the system at 
all "RG times" t is completely described by the single- 
conductance distribution function P(g\t). While being 
an uncontrolled approximation, EMA provides an instru- 
ment to attack the final stage of evolution of any initial 
distribution - the stage with a ~ g. We will see that, as 
it is typical for EMA, it works quite well, except for the 
immediate vicinity of MIT, for determination of energy- 
dependent effective conductance g g(t). 

Within the EMA, 



9ij 



(f-l)se ff 



(17) 



The effective conductance is then found from the self- 
consistency condition 



(Rij (gtj -5eff)) s 



0. 



(18) 



Thus, to find g c s(t) one should, in principle, solve 
Eqs. J3J) and (|17f) with an arbitrary given g c g(t) and 
find gijit) — g[gij{0), {gcfi}\t} as a functional of yet un- 
known function g c ff(t), and then, finally, solve Eq. 118|) 
for g c ff(t). This leads to a nonlinear integral equation 



p f s ■ 9[go, {9cs}\t] - 9cs{t) 
[go ) ago j Y\ +i i (z 7\ 777 



0. (19) 



For a general Po(ffo) this program can be fulfilled only 
numerically. If the distribution Pa{go) is narrow, an ex- 
plicit solution can be obtained for Sgij(t) — g%j{t) — ~g{t). 
For the standard deviation one finds ctema(s)/c(0) = 
e -s(l-2/z)_ c om p ar i son f this result (for the square lat- 
tice case z = 4) with the exact perturbation theory ljll|l 
is shown in Fig. At earlier stages (s <■ 1) agreement 
is rather good, but it becomes worse at large s where p- 
dependence of the eigenvalue A_ (p) becomes important. 




FIG. 4: Evolution of individual conductances according to 
the effective medium approximation. g c g(t) is shown with a 
thick line. 



Thus it seems that EMA may work reasonably good for 
broad distributions when s never becomes large. 

An important and physically relevant class of distribu- 
tions which allow for analytical EMA treatment is defined 
by the condition that In g is symmetrically distributed 
around some mean value. Writing the tunneling con- 
ductance as gij = ~g e~ hij one should require the distri- 
bution po(h) of fluctuations hij — n(dij — d), where dij 
are thicknesses of intcrgrain insulating barriers, to satisfy 
Po(h) — po(—h). For all such distributions on the square 
lattice (z = 4) a simple solution for the effective conduc- 
tance can be obtained: g e s(t) — g — t for t < t c = g 
and g c ff(t) = at t > t c . The individual conductances 
evolve as follows: for t < t c , 



gij(0)-g o , [^(o) + 5 ] 2 



(20) 



while for t > t c 



9v(t) = [ff«(* c ) - 2(t - t c )] %ij(t c ) - 2(t - t c )]. (21) 

This evolution is shown in Fig.^J The straight line g c g(t) 
is the separatrix: solutions g(t) with g(0) > ~g go above 
g e s(t) and eventually - one by one - become identical 
zeros at t > t c . Solutions with g(0) < ~g go below g e g(t) 
and become zeros all at once — at t = t c , together with 
g c s(t). For t c — t <C t c the latter solutions form a nar- 
row "bunch" , manifested by a sharp peak with the total 
intensity w 1/2 in the distribution function P(g,t) at 
g ~ (t c — t) 2 /g . For t > t c this narrow peak transforms 
into a J-peak at g = 0. Its intensity 1 — iV con( j(i) [with 
-^cond(i) being the fraction of conducting bonds, having 
gij > within the one-loop accuracy of Eq. (J2j] jumps 
from to 1/2 at t = t c and then grows monotonically, 
approaching unity at t 3> t c . 

The i-dependence of g e s(t) and N con( i(t) is shown in 
Fig. [S] together with the results of numerical simulations 
for the Cauchy initial distribution Po(g) — ^° 2 9 ^^ on 
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FIG. 5: Comparison of the EMA (dashed lines) and simu- 
lation (solid lines) results for the global conductivity g e s(t) 
and the fraction of conducting bonds iV con d(t) (see text for 
details) . 

the square lattice 32 x 32. One can see, that the simu- 
lated g c s{t) follows the EMA in the wide range of t, while 
in the vicinity of the transition it clearly deviates from 
the EMA and approaches zero with an exponent fi > 1. 
The numerically found N con ^(t) behaves, however, quite 
smoothly, showing no jump at t = t c . The reason for 
this smoothness is, apparently the tendency to clusteri- 
zation, demonstrated in Fig. ^ It can be interpreted as 
an instability, which drives especially small conductances 
(provided they form an appropriate configuration) to col- 
lapse earlier, than typical ones. 

Clustering and percolation. — Let us suppose that, due 
to a fluctuation in the initial distribution and/or to the 
dynamical evolution up to time t = to, the conduc- 
tances gij(to) of all bonds within a certain group hap- 
pen to be much smaller, than conductances in their im- 
mediate surrounding. The group should be such that 
the small conductances belong to the shell of some clus- 
ter of sites, separating it from the rest of the array (cf. 
Fig. Then, for any bond ij from the shell one gets 
R7 w Rg 1 = J2shci\9 kl - It means ! m particular, that 
the RG equations for the conductances of the shell are 
split from the rest of the RG equations, and can be solved 
separately. As a result, one obtains: 

,,s s (tshell ~ t) 0(*sheU - t) . 

9ij(t) ~ 9ij{to) , (22) 

tshell ^ £f) 

where i s hcii — to — -Re; 1 "C t c . Thus, we conclude, that, 
in contrast to the EMA solution, clusters surrounded by 
poorly conducting shells may become disconnected from 
the rest of the array already at t < t c . 

Numerical simulation clearly demonstrates the forma- 
tion of clusters (see Fig.0. When a disconnected cluster 
appears, the matrix A defined above Eq. acquires a 
new zero eigenvalue. Thus, the total number of electri- 
cally decoupled grains is given by the number of zero 
eigenvalues of the matrix A. The position t c of the 
percolative transition is a functional of the initial dis- 



tribution Po(g)- Apparently, t c is of the order of some 
mean initial conductance g , while the correct coefficient 
should be determined numerically. 

In general, evolution of the network of conductances 
with the growth of the parameter t{E) is rather similar 
to that would be expected at the classical percolation 
transition. However, our system can not be described by 
either purely "bond" or "site" percolation, due to devel- 
opment of local correlations (clustering) along with the 
RG flow. In particular, numerically observed (cf. Fig. 03) 
value of N con( i(t c ) is clearly larger than 1/2, contrary to 
expectations for purely bond percolation in 2D. More de- 
tailed numerical work is needed to determine the nature 
of this new kind of percolative transition; in particular, 
"measurements" of the conductivity exponent a (equal 
to 1.3 in the standard percolation problem |ljj) would 
be very desirable. 

Conclusions. — We have shown that at low temper- 
atures strong intrinsic inhomogeneities are developing 
in granular metal arrays with moderately large ran- 
dom bare conductances gij 3> 1. As a result, the 
Coulomb-driven metal-insulator transition expected if 
,9o < ln(Ec/5) 0, Q acquires features of percolation 
transition. Most directly the predicted behavior can be 
detected by measuring the distribution of the local tun- 
neling density of states at low temperatures. The best 
object for such a study would be a granular cermet of 
metal grains in the insulating matrix, like those studied 
in Refs. [LH Hj| . In these materials the ratio Ec/S was 
about 10 3 , indicating the existence of a broad range for 
logarithmic corrections to conductivity. It is hardly pos- 
sible that local tunnel conductances in such a granular 
cermet are all equal; at best, they can be distributed with 
the width of the order of the mean conductance. Our re- 
sults presented in Fig. [3] show that a simple logarithmic 
dependence g e s(T) — g — ln(g Ec/T) holds in a wide 
range of T for moderately random granular arrays as well, 
at least for the class of practically important symmetric 
distributions of log(<?) in the 2D space. 

If a granular metal has a tendency to become super- 
conductive with T sc ~ T c , its local superconductive prop- 
erties are expected to be strongly inhomogeneous due to 
position-dependent Coulomb effects. In other terms, su- 
perconductive properties of granular metal can be much 
more of "granular nature" than its normal properties at 
elevated temperatures. In this regard we mention very 
interesting recent experimental results |l6j |. 
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